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Preconditioning  techniques  that  are  used  to  alleviate  numerical  stiffness  due  to  low  Mach 
numbers  in  steady  flows  have  not  performed  well  in  unsteady  environments  since  the 
preconditioning  parameters  that  are  optimal  for  efficiency  are  detrimental  to  the  level  of 
spatial  dissipation  necessary  for  accuracy.  A  unified  flux  formulation  is  presented  where  the 
optimal  scaling  required  for  spatial  accuracy  is  independent  of  the  preconditioning  required 
for  convergence  thus  providing  a  framework  that  is  valid  over  a  broad  range  of  flow 
conditions.  Both  upwind  flux-difference  and  AUSM-type  schemes  are  investigated.  In  both 
cases,  the  use  of  unsteady  preconditioning  scaling  in  the  flux  formulation  is  shown  to  be 
critical  for  preserving  unsteady  accuracy.  In  the  flux-difference  case,  the  formulation  is 
based  on  a  generalized  blending  of  the  steady  and  unsteady  preconditioning  terms.  In  the 
AUSM  case,  the  formulation  introduces  two  modifications  to  the  standard  AUSM+m/? 
scheme,  designated  as  AUSM+m/U  wherein  the  pressure  dissipation  is  scaled  using  unsteady 
preconditioning  and  AUSM +u’p’  wherein  both  the  pressure  and  velocity  dissipation  terms 
are  scaled  by  the  unsteady  preconditioning.  Low  Mach  number  vortex  propagation  and 
acoustic  problems  are  used  to  demonstrate  the  strengths  of  the  formulation.  These  studies 
show  that  the  AUSM  family  generally  performs  better  than  the  blended  flux-difference 
schemes  in  terms  of  vortex  shape  preservation  and  control  of  odd-even  splitting  errors. 


I.  Introduction 

Accurate  and  efficient  modeling  of  unsteady  low  Mach  number  flows  continues  to  be  a  challenging  problem 
since  it  requires  the  resolution  of  disparate  time  scales.  Unsteady  effects  may  arise  from  a  combination  of 
hydrodynamic  effects  in  which  pressure  fluctuations  are  generated  by  unsteady  velocity  fluctuations,  and  acoustic 
effects  where  pressure  fluctuations  correlate  with  density  fluctuations  even  if  the  mean  Mach  number  is  very  low. 
Examples  in  the  former  category  would  include  vortex  propagation  problems  and  bluff  body  shedding,  while 
examples  of  the  latter  include  both  aeroacoustic  problems  as  well  as  acoustic  wave  propagation  in  internal  flows 
among  others.  Many  practical  applications  including  rotorcraft  flows,  jets  and  shear  layers  include  a  combination  of 
both  acoustic  and  hydrodynamic  effects.  Furthermore  these  effects  may  be  localized  with  the  flow  characteristics 
varying  in  the  domain.  Typically  the  near-field  may  exhibit  non-linear  coupling  between  the  different  modes  of 
unsteadiness,  while  the  far-field  may  show  a  separation  of  vortex  propagation  and  acoustic  scales  more  typical  of 
low  Mach  number  flows.  Therefore,  it  is  important  that  an  algorithm  designed  for  unsteady  low  Mach  number  flows 
function  efficiently  over  a  broad  range  of  flow  conditions  and  accurately  resolve  the  inherent  stiffness  of  the  system. 

Accurate  unsteady  simulations  of  these  complex  low  Mach  number  flows  pose  difficulties  for  time-marching 
algorithms  and  preconditioning  techniques  that  are  traditionally  used  to  alleviate  numerical  stiffness  work  well  for 
steady  state  simulations  but  have  problems  with  both  efficiency  and  accuracy  for  unsteady  computations.  One 
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crucial  factor  contributing  to  this  underperformance  arises  from  the  conflicting  requirements  placed  on  the  flux 
procedure:  preconditioning  parameters  that  are  optimal  for  time  scaling  (and  convergence)  are  detrimental  to  the 
level  of  spatial  dissipation  necessary  for  accurate  velocity  and  energy  transport.  The  goal  of  the  current  effort  is  the 
development  of  numerical  flux  procedures  where  the  optimal  scaling  required  for  spatial  accuracy  is  independent  of 
the  preconditioning  required  for  convergence  acceleration  thus  allowing  for  a  unified  formulation  that  would  be  both 
efficient  and  accurate  over  a  wide  range  of  flow  conditions. 

The  development  of  preconditioning  algorithms  has  been  a  very  active  area  of  research  in  the  literature  and  many 
flavors  of  flux  procedures  have  been  developed.  These  schemes  can  be  broadly  classified  into  two  families  based  on 
their  general  characteristics  for  low  Mach  number  flows:  1)  flux  difference/matrix  dissipation  preconditioning 
algorithms1"4  and  2)  the  AUSM  family  of  flux  split  schemes5'8.  The  flux-differenced  algorithms  have  been  applied 
with  much  success  to  steady  low  Mach  number  flows  by  recognizing  that  the  standard  flux  formulations  for  higher 
Mach  numbers  generate  too  little  dissipation  in  the  pressure  field  and  too  much  dissipation  in  the  velocity  field.  By 
multiplying  the  equations  by  a  preconditioning  matrix  that  effectively  scales  the  acoustic  speed  to  the  local 
convective  velocity,  optimal  spatial  dissipation  is  obtained  for  “all  speeds”  yielding  both  good  solution  convergence 
and  accuracy  for  steady  flows. 

For  unsteady  flows  characterized  by  high  local  cell  Strouhal  numbers  resolving  the  acoustic  time  scales  become 
important  and  the  steady  preconditioning  matrix  that  effectively  filters  the  acoustics  from  the  solution  becomes  far 
too  dissipative  for  resolving  the  pressure  wave  propagation.  If  more  generalized  definitions  of  the  preconditioning 
parameter  are  specified  that  take  this  local  Strouhal  number  into  account9,  the  preconditioning  parameter  reverts 
back  to  its  unpreconditioned  value  as  acoustic  effects  become  dominant  thus  restoring  the  ability  to  accurately 
model  acoustic  propagation  with  good  convergence.  However,  the  key  drawback  of  flux  differenced  procedures  is 
that  since  the  preconditioning  parameter  used  for  time  scaling  also  affects  the  dissipation  for  the  spatial  flux, 
hydrodynamic  unsteady  effects  (such  as  vortex  propagation)  where  the  velocity  and  temperature  fields  are  still 
governed  by  the  convective  time  scales  become  inaccurate  as  the  preconditioning  parameter  reverts  to  its  original 
acoustically  accurate  form.  Therefore,  solutions  to  low  Mach  number  vortex  propagation  problems  may  show  a 
disconcerting  disconnect  between  convergence  and  accuracy  of  the  velocity  field.1 

In  an  effort  to  improve  on  this  behavior,  Sankaran  and  Merkle  reformulated  the  flux  difference  scheme  with 
features  akin  to  those  of  AUSM  schemes.9  Here,  the  eigenvalue  matrix  used  for  spatial  flux  dissipation  was  altered 
such  that  the  eigenvalue  associated  with  pressure  equation  corresponded  to  the  acoustic  wave  from  the 
preconditioned  system  whereas  all  the  other  equations  use  the  convective  velocity.  While  accurate  results  with 
unsteady  preconditioning  are  reported  for  the  vortex  propagation  problem  the  concern  with  this  formulation  is  that 
the  spatial  dissipation  does  not  reduce  to  the  “steady”  preconditioning  form  if  unsteady  effects  are  not  dominant.  To 
avoid  this  discrepancy  for  steady  calculations,  Potsdam  et  al1.  used  a  “blended”  procedure  by  defining  selection 
matrices  that  use  the  flux  from  the  “unsteady”  preconditioning  matrix  for  the  continuity  equation  and  the  flux  from 
“steady”  preconditioning  for  the  momentum  and  energy  equations.1  In  the  limit  where  the  flow  is  steady  the 
“unsteady”  preconditioning  parameter  for  the  continuity  equation  reverts  back  to  the  steady  preconditioning 
parameter  and  the  overall  steady  preconditioning  formulation  is  unchanged. 

The  second  class  of  numerical  flux  procedures,  the  AUSM  family  of  flux-split  schemes  has  also  been  extended 
5-8 

to  low  Mach  number  flows  .  A  key  characteristic  of  AUSM  schemes  is  that  the  convective  flux  is  separated  from 
the  pressure  flux;  both  flux  terms  are  computed  independently  as  scalar  formulations  thus  making  it  possible  to 
independently  tailor  dissipation  for  hydrodynamic  and  acoustic  unsteadiness.  The  convective  flux  terms  have 
essentially  a  upwinded  form  with  additional  dissipation  at  low  Mach  numbers;  Liou  demonstrates  excellent  results 
for  steady  low  Mach  number  flows  in  the  AUSM+i/p  formulation.5  The  pressure  flux  term  also  has  additional 
dissipation  terms  in  the  AUSM+i/p  formulation  but  these  are  aimed  at  resolving  instabilities  at  the  sonic  point.  In 
addition  to  the  AUSM+wp  formulation  we  will  also  focus  on  a  variant  of  the  AUSM  scheme  by  Shima  and 
Kitamura7,8  which  is  referred  to  as  the  SLAU  scheme  (Simple  Low  Dissipation  AUSM)  by  the  authors  and  has  been 
applied  primarily  to  unsteady  flows.  A  detailed  discussion  comparing  the  AUSM+wp  and  SLAU  scheme  will  be 
presented  in  the  paper  and  the  SLAU  scheme  is  shown  to  be  a  equivalent  to  the  original  AUSM+  scheme  in  the  low 
Mach  number  limit,  i.e.,  without  inclusion  of  the  pressure  and  velocity  dissipation  terms. 

The  primary  focus  of  this  paper  is  to  present  a  more  generalized  preconditioning  framework  based  on  an 
unsteady  Mach  number  parameter  which  ensures  that  the  flux  formulation  is  accurate  and  efficient  for  both 
hydrodynamic  and  acoustic  instabilities,  and  reverts  to  traditional  steady  flux  form  in  the  limit  of  steady  low  Mach 
number  flows.  The  generalized  unsteady  preconditioning  framework  has  been  adapted  for  both  AUSM  and  flux- 
differenced  family  of  schemes.  New  flux  differenced  procedures  have  been  derived  by  developing  “blended” 
formulations  that  mix  “unsteady”  and  “steady”  preconditioning  parameters  locally  for  each  term  in  the  dissipation 
matrix  and  are  extensions  of  the  work  presented  by  Potsdam  et  al1.  For  the  AUSM  family  of  schemes  the 
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dissipation  for  both  the  mass  flux  and  pressure  terms  have  been  modified  using  the  unsteady  Mach  number 
parameter.  The  modified  AUSM  schemes  are  referred  to  as  AUSM+i/p'  and  AUSM+w/U,  where  the  “primes” 
designate  that  the  pressure  and/or  velocity  dissipation  terms  are  being  scaled  by  the  unsteady  Mach  number. 
Corresponding  modifications  are  proposed  for  the  SLAU  schemes  as  well.  Both  the  blended  flux-difference  and  the 
modified  AUSM+wp  formulations  are  tested  on  a  range  of  unit  problems  encompassing  both  hydrodynamic  and 
acoustic  instabilities.  Both  solution  accuracy  and  unsteady  inner  iteration  convergence  are  evaluated  for  these 
problems  to  demonstrate  the  wide  range  of  applicability. 

II.  Flux  Formulations  for  Low  Mach  Number  Unsteady  Flows 

The  formulation  of  a  more  generalized  unsteady  preconditioning  framework  for  both  flux  differenced  and  AUSM 
family  of  schemes  is  described  here.  We  begin  by  giving  a  brief  overview  of  the  preconditioned  system  of  equations 
and  steady  preconditioning  for-flux  difference  schemes  in  Sections  II. A  and  II.B.  The  “blended”  flux  difference 
formulations  for  unsteady  preconditioning  that  were  developed  as  part  of  this  effort  are  described  in  Section  II. C. 
The  AUSM  family  of  schemes  are  described  in  Section  II. D  and  an  analysis  comparing  AUSM+up  and  SLAU  is 
discussed  here.  The  extensions  for  AUSM+i/p '  and  AUSM+i/  ’p  ’  employing  the  unsteady  preconditioning  parameter 
are  described  in  Section  II. E. 


A.  Preconditioned  Equation  System 

The  standard  conservative  form  of  the  equations  (with  a  two -equation  turbulence  model)  may  be  written  as 

dQ  dF  dF  dG  C|T^ 

—  +  —  +  —  + — =S+DV 
dt  dx  dy  dz 


(1) 


where  Q  =  [p,  pu,  pv,  pw,  e,  pkps]1  represents  the  vector  of  conserved  variables;  E ,  F  and  G  are  the  flux  vectors;  S 
represents  the  source  terms  (if  any);  and  Dv  represents  the  viscous  fluxes.  For  low  Mach  number  flows  as  density 
approaches  a  constant  the  conservative  variables  become  ineffective  for  temporal  integration.  In  this  flow  regime, 
the  primitive  variables  vector,  Qv  =  [p,  u,  v,  w,  T,  k,  s]J ,  constitutes  a  particularly  effective  choice.  Primitive 
variables  replace  the  density  by  pressure,  thereby  avoiding  round-off  errors,  and  the  total  energy  by  temperature 
which  is  more  compatible  with  the  thermal  diffusion  terms  and  for  modeling  multi-species  flows  (as,  for  example, 
the  combusted  exhaust  plume  from  an  aircraft  engine). 

One  effective  way  for  expressing  a  general  iterative  method  is  through  a  dual -time  formulation.  Upon  appending 
a  pseudo-time  derivative  and  a  preconditioning  matrix,  Tp,  Eqn.  (1)  takes  the  form; 


-j— ,  dQv  dQ  dE  dF  dG 

r-^+-^ + — + — + —  =  s+dv 

dr  dt  dx  dy  dz 


(2) 


The  numerical  characteristics  of  the  pseudo-time  in  this  equation  are  determined  by  the  eigenvalues  of  the  matrix, 
[r;1  {dE!dQvj\  which  are: 

~1 


A  =  diag 


2L 


U(l  +  R)±ylu2(  1-tf)2  +4c'2(/;  +/’  +/z2) 


,u,u,u,u,u 


(3) 


where  R  represents  the  ration  of  the  square  of  the  pseudo-acoustic  speed  and  the  physical  acoustic  speed  which  can 
be  expressed  as  the  square  of  the  ratio  of  the  physical  Mach  number  to  a  pseudo-Mach  number  as 


R  = 


(ppPhT  +  Pj (i  pfop)} 


M_ 

aF 


(4) 


{ppPhT  Pt pfopl^j  v  £  j 

This  system  of  equations  can  be  made  well-conditioned  by  choosing  pp  to  be  of  order  p'p  «^0(l/w2)J,  thereby 

resulting  in  a  pseudo-  sound  speed,  d  « [0(w)] .  The  preconditioning  parameter,  pp  ,  can  be  expressed  in  terms  of 
the  ratio  of  pseudo-Mach  number  to  physical-Mach  number  and  the  ratio  of  specific  heats: 

\2 


Pp 

Pp 


I  M' 

V  \mp 
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where 


M]  =min[max(V2,V/V;,n),l] 


(5) 


(6) 


3 

American  Institute  of  Aeronautics  and  Astronautics 


Here,  Mt  is  the  local  Mach  number,  Mmin  is  a  user  specified  cut-off  Mach  number  to  preclude  difficulties  at 
stagnation  points,  and  Mu  is  an  “unsteady”  preconditioning  scaling.  This  unsteady  preconditioning  term  can  be 
related  to  the  global  Strouhal  number  for  the  problem  as: 


Mu  = 


L 

nc/St 


(7) 


where  L  is  a  global  length  scale  and  At  is  the  time-step. 


B.  Flux-Difference  Algorithm:  Roe’s  Scheme  with  Steady  or  Unsteady  Preconditioning 

The  inviscid  flux  at  cell  interfaces  is  calculated  with  an  approximate  Riemann  solver  given  the  reconstructed  left 
and  right  solution  states.  The  flux  difference/matrix-dissipation  procedure  based  on  Roe’s  scheme11  can  be  given  as: 

Fm  =  ~  |  F(Qv ,  n) + F(Qy ,  n) + AF”  (Qv ,  n) + AF+  (Qv ,  n)  j  (8) 

where  n  is  the  area-directed  normal  of  either  the  cell  face  or  the  dual  face  crossing  the  edge,  Q“  and  are  the  left 
and  right  primitive  state  variables,  and  Qv  denotes  the  Roe-averaged  variables.  The  flux  differences  AF"  and  AF+ 
are  given  by: 

AF+(Qv,noi)  =  i (fpRp(Ap+|  Ap  |)Lp)(Q+-Q-)  (9) 

AF- (Qv , noi )  =  -f  (f pRp (Ap - 1  Ap  |)Lp)(Q+-Q-)  (10) 

Here  rp,  Rp,  Lp,  and  Ap  are  the  preconditioning  matrix,  right  eigenvector  matrix,  left  eigenvector  matrix,  and 
diagonal  eigenvalue  matrix  for  the  preconditioned  system  computed  with  the  Roe-averaged  variables. 

C.  Flux-Difference  Algorithms:  Blended  Unsteady/Steady  Preconditioned  Schemes 

The  standard  flux  difference  formulation  shown  above  (Eqns.  (9)  and  (10))  indicate  that  the  spatial  dissipation  is 
intimately  tied  to  the  time  scaling  defined  by  the  preconditioning  matrix.  Therefore,  if  “unsteady”  preconditioning 
is  used,  as  given  by  Eqn.  (7),  for  unsteady  low  Mach  number  flows  good  inner  iteration  convergence  is  obtained  at 
the  expense  of  solution  accuracy  and  is  clearly  not  optimal.  To  develop  flux  differenced  procedures  that  can  provide 
accurate  unsteady  solutions  while  continuing  to  provide  good  convergence,  generalized  “blending”  formulations 
were  constructed  in  which  “steady”  and  “unsteady”  preconditioning  forms  could  be  used  for  specific  terms  in  the 
momentum  and  energy  equations.  These  formulations  are  devised  by  exploiting  the  form  of  the  algebraic 
expressions  for  the  spatial  dissipation  terms.  The  combination  of  the  spatial  dissipation  terms  in  Eqns.  (9)  and  (10) 
may  be  written  as  a  matrix,  Cp,  multiplied  by  the  change  in  the  primary  dependent  variable  vector: 

A^(Qv,fioi)+AF-(Qv,noi)  =  Cp  (<#-<£)  0D 

where  the  combined  dissipation  matrix,  Cp,  can  be  written  out  as 


C11 

C12 

Cl3 

CM 

Cl5 

Cie 

C17 

uCll+C21 
yC|  i  +C3] 

uC12+C22 

uC13+C23 

uC14+C24 

uC15  +C25 

uC16 

uCi? 

(12) 

wC'n  +C41 
HCn+Cjj 

hc12+<4 

HC13+C53 

HC14  +C54 

hc15+c;5 

hc16 

hc17 

kCn 

kC12 

kC13 

kC14 

kC15 

kC16 

kCi7 

sCu 

... 

"  *  ) 

The  form  of  this  matrix  shows  that  the  entries  for  the  scalar  transport  equations,  such  as  those  for  the  turbulence 
scalars  or  species  equations  for  combusting  flows,  are  multiples  of  the  first  row  (the  continuity  equation)  scaled  by 
the  convected  property  of  that  equation.  Similarly,  the  entries  for  the  momentum  and  energy  equations  contained 
one  term  that  is  a  product  of  the  first  row  multiplied  by  the  pertinent  local  scalar  (e.g.,  u,  H),  however,  they  also 
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include  an  additive  correction  term  which  is  a  function  of  the  eigenvalues.  Therefore,  the  form  of  this  matrix 
suggests  two  “blended”  schemes  in  which  steady  and  unsteady  preconditioning  can  be  applied  to  different  parts  of 
this  matrix. 

The  first  blended  scheme  was  constructed  by  decomposing  the  dissipation  matrix  into  two  terms  as  given  by 

AF++AF-=(cOT+C^)(q+-Q-)  (13) 


where  the  first  term  is  given  by  the  tensor  product  of  two  vectors  given  by  C-  (Cn  C12  C13  C14  C15  C16  C17 ) 

and  0  =  (l  u  v  w  H  k  s)  and  the  second  term  is  a  correction  matrix  C' .  This  matrix  structure  suggests  a 

sophisticated  blending  strategy  in  which  “steady”  and  “unsteady”  preconditioning  forms  could  be  used  for  specific 
terms  in  the  momentum  and  energy  equations  which  may  be  written  as: 

AF+  +  AF-  =|c<D  +C^)(Q+-Qv)  (14) 

Here  the  COr  terms  use  “unsteady”  (green  notation)  while  the  Cp  terms  could  use  “steady”  preconditioning  (red 

notation).  This  blending  scheme  is  equivalent  to  having  the  first  row  of  the  Cp  matrix  with  unsteady  preconditioning 
while  the  other  rows  will  have  a  combination  of  steady  and  unsteady  preconditioning  terms. 

Another  blended  formulation  can  be  constructed  in  which  the  “unsteady”  eigenvalue  and  preconditioning  matrix 
is  used  to  evaluate  the  pressure  equation  while  the  corresponding  “steady”  values  are  used  to  evaluate  the  flux  for 
the  momentum  and  energy  equations: 


AF+(Qv,iioi)  = 


lp  + 


( 


rp,sRp,s(A+lAl)p,sL. 


p,s)Lv,t] 


(Ql-QJ 


(15) 


This  blending  scheme  is  equivalent  to  having  the  first  column  of  the  Cp  matrix  in  Eqn.  (12)  with  “unsteady”  (green 
notation)  preconditioning  terms  while  the  other  columns  will  be  evaluated  with  “steady”  (red  notation) 
preconditioning  terms.  The  selection  matrices  in  Eqn.  (15),  Lp  and  LV?T,  are  used  to  select  the  continuity  equation  or 
momentum  and  energy  equations.  This  blended  scheme  is  very  similar  the  scheme  proposed  by  Potsdam  et  al. 
which  is  given  by: 1 


^  (Qv  ’  aoi )  =  \  [(rp  A.»  (A+  I  A  |)p>u  Lp  u  )  Lp  + 
(fpA,s(A+lAl)p,sLp,s)Lv,T 


(q;-q;) 


(16) 


The  sole  difference  between  the  two  schemes  is  that  the  preconditioning  matrix  in  the  second  set  of  terms  in  Eqn. 
(15)  uses  steady  preconditioning  (red  color)  while  Potsdam’s  method,  Eqn.  (16),  uses  unsteady  preconditioning 
(green  color).  However,  the  form  of  this  flux  formulation  cannot  be  reduced  to  altering  the  first  column  of  the  Cp 
matrix  as  was  done  for  the  blended  method  given  by  Eqn.  (15). 

In  summary,  using  the  exact  algebraic  form  of  the  flux-difference  spatial  dissipation  term  allowed  for  the 
development  of  two  generalized  blended  flux  formulations  which  feature  a  combination  of  unsteady  preconditioning 
for  the  continuity  equation  and  steady  preconditioning  for  the  momentum  and  energy  equations.  However,  the 
solutions  for  various  unsteady  low  Mach  number  test  cases  predicted  by  these  formulations  as  well  as  with 
Potsdam’s  blended  flux  scheme  were  indistinguishable.  Therefore,  numerical  results  for  only  one  of  these  schemes 
will  be  presented  in  the  Section  III. 


D.  AUSM  Type  Algorithms  for  Preconditioned  Systems 

Flux-splitting  schemes  provide  an  alternative  approach  to  flux  difference  schemes  for  calculating  the  inviscid 
flux  at  cell  interfaces.  These  procedures  are  attractive  since  they  provide  a  natural  means  to  decouple  the  dissipation 
for  the  momentum  and  pressure  flux  thus  allowing  more  flexibility  to  tailor  the  dissipation  for  acoustic  and 
hydrodynamic  instabilities  in  unsteady  flows.  Prior  to  describing  the  extensions  to  an  unsteady  preconditioning 
framework  we  analyze  the  differences  between  the  AUSM+i/p  scheme5  developed  by  Liou  and  the  SLAU  scheme 
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derived  from  AUSM  by  Shima  and  Kitamura7,8  specifically  for  low  Mach  number  flows.  These  schemes  are 
described  below. 

The  numerical  flux  for  flux-split  schemes  is  given  by 

~  m  +  \m\  ^  m-\m 

2  2 

0-(l,u,v,w,hr)T  (18) 


-0  +pn 


(17) 


The  unique  formulation  of  the  AUSM  family  schemes  is  determined  by  the  choice  of  mass  flux  function,  ill ,  and 
average  pressure,  p  .  The  mass  flux  in  AUSM +up  scheme  is  given  by 


m  =  cMV2 


\pL  f  M\n  >  0 
[pR  otherwise  ’ 


MV2  =  )  +  M~4)  (Ms  )  -  ^  rnax(l  -  aM: 

J  a 


,0)^ 

pc 


(19) 


Here  are  fourth  order  polynomials  in  Mach  number.  The  last  term  is  a  dissipation  term  formulated  for  low 

Mach  number  flows.  The  fa  term  in  the  denominator  increases  the  dissipation  in  inverse  proportion  to  the  Mach 
number  and  is  equivalent  to  a  steady  preconditioning  parameter.  It  is  defined  as: 

fa  =  M0 (2 — M0 ) e [0, 1],  /W2=min(l,max(V2,^lln)),  M2  =(u2L+u2)/ 2c2  (20) 

Making  a  low  Mach  number  assumption  and  dropping  higher  order  terms  in  Mach  number  Eqn.  (19)  can  be 
rewritten  as: 

m  =  /r  +  ^")-^max(l-<rM2, 0)^33^!  (21) 

where  p±  is  p+  for  Mx/ 2  greater  than  zero  and  p  otherwise. 

The  pressure  flux  formulation  in  the  AUSM+i/p  scheme  is  given  by; 

P  =  PsfML)PL  +]T5)(Mr)Pr  ~2KuPP5pffSLl  (22) 


Here  are  fifth-order  polynomials  in  Mach  number  while  the  last  term,  which  can  be  summarized  as  the  pu 

dissipation  term,  is  a  dissipation  term  that  operates  at  transonic  Mach  numbers  and  decreases  in  value  as  the  Mach 
number  drops.  Making  a  low  Mach  number  assumption  and  dropping  higher  order  terms  in  Mach  number  this 
equation  reduces  to: 


=  ^  (p+  +  P~  )  ^ 1  (M+P+  -  M~p~  )-lv„[2M](2- 3M~  +  3 M+  )  pcAu  (23) 


where  the  [2 M\  factor  in  the  last  term  comes  from  the  expression  for^. 

For  the  SLAU  scheme  given  by  Shima  and  Kitamura7,8,  the  mass  flux  is  given  directly  in  the  simplified  low 
Mach  number  form  as 


*=&(*? +\K\yp-{K-p'  r)}-#4? 

%  =  1-2M 

While  the  pressure  flux  for  the  SLAU  formulation  is  given  by 


(p++p )  + 


(5)  (5) 


(P+  -  P~)  +  (1 ' -  Z)(P5)  +p5)~l) 


C P++P~) 


(24) 


(25) 


Making  a  low  Mach  number  assumption  and  dropping  higher  order  terms  in  Mach  number  this  equation  reduces  to: 


+p-)+^(M+  +M-)(p+  +p)~l[2M]pcAu 


(26) 
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where  the  [2 M]  factor  in  the  last  term  comes  from  the  expression  for  %. 

Comparing  the  mass  flux  dissipation  for  AUSM+i/p  (Eqn.  (21))  and  SLAU  (Eqn.  (24)),  it  is  observed  that  the 
primary  difference  for  the  two  schemes  is  the  dissipation  that  is  added  at  low  Mach  numbers  for  the  AUSM+i/p 
scheme.  In  the  AUSM+i/p  scheme  the  dissipation  term  for  the  mass  flux  (which  is  a  pressure  dissipation)  has  the^ 
term  in  the  denominator  which  increases  the  dissipation  substantially  as  the  Mach  number  drops  and  behaves  like  a 
steady  preconditioning  parameter.  In  contrast,  the  SLAU  formulation  does  not  have  this  term  and  exhibits  only  a 
weak  dependence  on  Mach  number  from  the  %  term  in  the  numerator.  Therefore,  the  SLAU  scheme  would  be 
inadequate  for  steady  low  Mach  number  problems  and  this  is  confirmed  by  demonstrating  the  formulation  on  a 
steady  test  case  in  Section  III  A. 

In  contrast,  the  pressure  flux  dissipation  (which  is  a  velocity  dissipation)  between  the  AUSM+wp  (Eqn.  (23))  and 
SLAU  (Eqn.  (25))  are  not  very  different  at  low  Mach  numbers.  The  dissipation  term  pu  in  the  AUSM+up 
formulation  is  designed  to  operate  close  to  the  sonic  point  and  drops  in  proportion  to  the  Mach  number.  However, 
as  will  be  described  later  for  very  low  Mach  number  acoustic  unsteady  problems,  additional  dissipation  is  needed  for 
the  pressure  flux  as  well.  This  dissipation  term  is  described  as  part  of  the  unsteady  extensions  in  the  following 
section  and  demonstrated  for  a  test  problem  in  Section  III. 


E.  Extensions  of  Unsteady  Preconditioning  Framework  to  AUSM  schemes 

As  our  analysis  in  the  previous  section  indicates  the  standard  AUSM +up  scheme  is  optimized  for  steady  low 
Mach  number  problems.  A  more  general  unified  formulation  that  automatically  selects  the  appropriate  dissipation 
level  is  presented  here  by  formulating  the  mass  flux  dissipation  in  terms  of  an  unsteady  Mach  number  scale.  This 
modified  scheme  is  referred  to  AUSM+wp  where  the  pressure  prime  indicates  that  the  pressure  dissipation  is  scaled 
by  the  unsteady  Mach  number.  Furthermore  the  pressure  dissipation  term  is  also  modified  by  adding  additional 
dissipation  for  low  Mach  number  acoustic  problems  using  the  unsteady  Mach  number  parameter.  These 
modifications  lead  to  the  AUSM+w ’p*  scheme  where  the  primes  indicate  that  both  pressure  and  velocity  dissipation 
terms  are  scaled  by  the  Mach  number  parameter. 

The  dissipation  for  the  mass  flux  term  in  the  AUSM+w/U  scheme  is  generalized  by  modifying  the  definition  of 
Mach  number  used  to  define  the  fa  term  in  Eqn.  (20)  by  including  an  unsteady  preconditioning  scale  as  follows: 

M02=min(l,max(M2,ML,M„2))  (27) 


We  note  that  in  the  most  general  configuration  the  Mu  scale  can  be  defined  both  as  a  global  parameter  and  as  local 
factor  that  is  computed  on  the  local  flow  physics.  For  simpler  problems  a  global  parameter  might  suffice  but  for 
more  complex  physics  where  the  characteristics  may  vary  (e.g.  high  frequency  in  the  near  field  of  a  turbulent  jet  and 
low  frequency  in  the  far-field)  the  ability  for  the  numerical  formulation  to  select  the  automatically  select  the 
unsteady  scale  is  crucial.  This  issue  is  discussed  further  in  the  following  section.  We  note  that  reducing  the 
dissipation  for  the  mass  flux  term  for  unsteady  acoustic  flows  has  been  suggested  both  by  Liou5and  Vigneron  et  al6 
however  Eqn.  (27)  provides  a  means  for  the  numerical  formulation  to  self-select  the  appropriate  dissipation  level. 

The  dissipation  for  the  pressure  flux  term  for  the  AUSM+ipU  scheme  has  also  been  modified  for  unsteady  low 
Mach  numbers.  The  definition  of  Mach  number  M0  using  the  unsteady  Mach  number  scale  Mu  as  shown  above  in 
Eqn.  (27)  also  affects  the  pressure  dissipation  term  in  Eqn.  (22)  as  shown  below: 


fa  &  0  when  Mu  for  steady  low  Mach  number  flows 
fa  ~  1  when  Mu&  \  for  unsteady  low  Mach  number  flows 


(28) 


In  particular,  for  steady  flows  no  additional  dissipation  results  but  for  unsteady  acoustic  low  Mach  number  flows 
substantial  dissipation  is  added  to  the  pressure  flux.  The  effect  of  this  additional  dissipation  for  acoustic,  low  Mach 
number  flows  will  be  demonstrated  on  a  test  problem  in  Section  HID. 

We  further  note  that  the  SLAU  scheme  in  fact  corresponds  to  the  unsteady  low  Mach  limit  of  AUSM+i/p  \  The 
unsteady  Mach  number  parameter  in  the  latter  formulation  naturally  decreases  the  pressure  dissipation  in  the 
acoustic  limit,  which  leads  to  the  “non-preconditioned”  formulation  that  is  inherent  in  SLAU.  Thus,  as  formulated, 
the  SLAU  scheme  appears  more  appropriate  for  unsteady  problems,  although  the  proper  pressure  flux  dissipation 
can  be  introduced  to  extend  it  to  steady  low  Mach  number  flows. 


III.  Evaluation  of  the  Low  Mach  Number  Formulations 

The  low  Mach  number  formulations  described  in  the  previous  section  were  evaluated  with  several  test 
problems  involving  hydrodynamic  and  acoustic  unsteadiness.  These  calculations  were  performed  with  the 
CRUNCH  CFD®  code,  developed  at  CRAFT  Tech12"15.  The  candidate  flux  formulations  for  unsteady  low  Mach 
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number  flows  will  be  tested  out  rigorously  for  the  following  three  test  cases  that  include  both  hydrodynamic  and 
acoustic  instabilities:  1)  Unsteady  inviscid  Lamb  vortex  problem  (hydrodynamic  instability),  2)  Unsteady  inviscid 
flow  in  a  pipe  with  fluctuating  back  pressure  (mixed  acoustic  and  hydrodynamic  instability),  and  3)  Shock  tube  with 
small  pressure  difference  (pure  acoustic  problem).  We  note  that  both  the  “blended”  flux  difference  and  AUSM  type 
schemes  have  to  be  run  with  inconsistent  LHS  and  RHS  discretization  since  these  formulations  are  unstable  when  a 
consistent  LHS  is  used.  For  the  results  presented  here  we  employ  the  standard  preconditioned  Roe  flux  differenced 
procedure  on  the  LHS  with  unsteady  preconditioning  and  dual  time  stepping  at  each  time  step.  Prior  to  discussing 
the  unsteady  test  cases,  a  steady  low  Mach  number  test  case  is  presented  to  compare  the  Roe  flux  differenced 
procedure  with  the  SLAU  scheme;  as  discussed  earlier  the  dissipation  in  the  original  SLAU  scheme  is  shown  to  be 
appropriate  only  for  unsteady  flows  and  this  constraint  is  demonstrated  by  computing  a  steady  flow- field. 


(a)  Pressure  contours  with  SLAU  Scheme:  Odd-even  oscillations  are  evident 


(b)  Pressure  coefficient  (c)  Residual  convergence 

Figure  1.  Solution  for  a  steady  flow  over  a  NACA0015  airfoil  at  a  freestream  Mach  number  of  0.001  and  an 
angle  of  attach  of  4°  as  computed  with  the  SLAU  scheme  and  a  Roe  flux  differenced  scheme  with  steady 
preconditioning. 

A.  Steady  Flow  over  a  NACA0015  Airfoil 

Steady  state  calculations  were  performed  for  a  NACA0015  airfoil,  at  a  freestream  Mach  number  of  0.001, 
Reynolds  number  of  1.95E+06,  and  an  angle-of-attack  of  4°  using  a  Roe  flux  differences  procedure  with  steady 
preconditioning  and  the  SLAU  scheme.  The  higher  order  solution  and  convergence  is  shown  Figure  1  for  both 
schemes.  For  the  SLAU  scheme  it  is  noted  that  the  solution  shows  oscillations  with  odd-even  coupling  that  is 
reflected  in  the  surface  pressure  coefficient.  Moreover,  the  corresponding  convergence  shown  in  Figure  1(c) 
indicates  that  the  SLAU  solution  does  not  converge  in  contrast  to  the  solution  with  flux  difference  procedure  with 
steady  preconditioning.  The  numerical  instability  for  the  SLAU  was  verified  to  result  from  the  low  freestream  Mach 
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numbers;  the  odd-even  oscillations  were  not  evident  when  the  Mach  number  was  increased  to  0.1  (results  not 
shown).  Therefore,  it  was  concluded  that  the  mass  flux  dissipation  in  Eqn.  (24)  for  SLAU  scheme  was  unable  to 
provide  accurate  solutions  for  steady  low  Mach  number  flows  as  expected. 

The  dissipation  term  for  SLAU  was  subsequently  modified  by  adding  the  form  of  the  dissipation  from  the 
AUSM+wp  formulation  (Eqn.  (21)).  Therefore,  for  steady  low  Mach  number  problems  substantial  dissipation 
(resulting  from  thefa  term  in  the  denominator  of  Eqn.  (21)  is  added  to  the  mass  flux  term,  leading  to  what  we  refer 
to  as  the  SLAU+p  scheme.  The  solution  for  the  NACA0015  airfoil  with  the  modified  dissipation  term  is  shown  in 
Figure  2.  The  flow  contours  and  pressure  coefficient  are  now  smooth  and  the  convergence  was  found  to  be 
independent  of  the  Mach  number.  Therefore,  this  confirms  our  premise  discussed  earlier  that  the  SLAU  scheme 
corresponds  to  the  standard  AUSM+  formulation  and  requires  modifications  in  the  steady  limit  (i.e.,  the  SLAU+p 
scheme)  for  low  Mach  number  flows. 


(a)  Pressure  contours  with  additional  dissipation  at  low  Mach  number  number  (SLAU+p  scheme) 


(b)  Pressure  coefficient  (c)  Residual  convergence 

Figure  2.  Solution  for  a  steady  flow  over  a  NACA0015  airfoil  at  a  freestream  Mach  number  of  0.001  and  an 
angle  of  attach  of  4°  as  computed  with  the  SLAU+p  scheme. 


B.  Vortex  Propagation  Problem 

The  propagation  of  an  unsteady  convecting  inviscid  Lamb  vortex  was  used  to  assess  the  effectiveness  of  the  flux 
formulations  discussed  above  for  the  hydrodynamic  class  of  problems.  The  velocity  distribution  in  polar 
coordinates  for  the  Lamb  vortex  is  given  by 


Vr=0 


ve  =  r 


\-e-r2/ 


(29) 


V 


y 
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where  F  and  §  are  the  vortex  strength  and  characteristic  radius  of  the  vortex  and  were  set  to  0.1  and  0.03, 
respectively.  Simulations  are  presented  here  for  a  freestream  Mach  number  of  0.001  were  computed  on  a  uniform 
Cartesian  grid  with  a  spacing  of  0.005.  The  vortex  solution  is  presented  after  the  vortex  has  travelled  a  distance 
equal  to  0.4  on  this  uniform  mesh.  The  unsteady  Mach  number  scale  was  specified  as  done  by  Potsdam  et  al.  in  Ref. 
Ref.  1. 


M=- 


L(MJ 


7r(Ax)(CFLu ) 


(30) 


For  a  baseline  value  of  CFLU  =  1  this  translates  to  an  unsteady  Mach  number  of  2.55E-02.  For  very  small  time  steps 
(e..g  CFLU  =  0.001)  the  unsteady  Mach  number  becomes  1  and  the  dissipation  levels  are  representative  of  unsteady 
preconditioning.  Results  are  presented  for  the  baseline  time-step  dictated  by  CFLU  =  1  and  an  extremely  small  time- 
step  given  by  CFLU  =  0.001.  Two-hundred  and  fifty  inner  iterations  were  used  for  all  cases.  Both  the  blended  flux 
difference  scheme  and  the  AUSM+wp'  scheme  are  compared  with  benchmark  values  of  Roe  flux  differenced 
scheme  using  steady  and  unsteady  preconditioning. 

The  vortex  solution  and  inner  iteration  convergence  are  shown  for  the  baseline  time  step  of  CFLU  =  1  in  Figure  3 
for  the  different  schemes.  The  exact  vorticity  profile  is  given  by  the  black  line  in  Figure  3(a).  The  red  and  blue 
lines  correspond  to  the  standard  flux-difference  scheme  with  steady  and  unsteady  preconditioning,  respectively.  The 
result  computed  with  a  blended  flux  formulation  is  given  by  the  green  line  whereas  the  orange  line  gives  the  profile 
computed  with  the  AUSM+i/p'  scheme.  The  vorticity  profiles  show  that  the  flux  difference  scheme  with  steady 
preconditioning,  the  blended  flux  difference  scheme,  and  the  AUSM+i/p'  collapse  to  essentially  the  same  solution. 
The  solution  computed  using  the  flux  difference  scheme  with  unsteady  preconditioning  has  poor  accuracy  compared 
to  the  other  schemes,  however,  the  convergence  plots  shows  that  it  has  the  best  rate  of  convergence  indicating 
optimal  time  scaling  for  convergence.  Rectification  of  this  inconsistency  between  solution  accuracy  and 
convergence  was  one  of  the  primary  goals  at  the  outset  of  this  research.  The  convergence  for  the  blended  scheme  is 
slightly  better  than  the  flux  difference  scheme  with  steady  preconditioning.  The  convergence  for  the  AUSM+wp' 
scheme  is  worse  than  the  blended  scheme  which  is  likely  due  to  the  increased  inconsistency  of  the  left-hand  side  and 
right-hand  side  flux  formulations,  however,  the  solution  was  just  as  accurate  as  the  blended  and  pure  steady 
preconditioning  schemes. 


(b)  Inner  iteration  convergence 


Figure  3.  Vortex  solution  accuracy  and  convergence  a  time-step  of  CFLU  =  1. 


To  further  confirm  the  insensitivity  of  the  solution  accuracy  to  the  unsteady  preconditioning  parameter,  the  time- 
step  was  reduced  to  CFLU  =  0.001  which  goes  to  the  limit  of  unsteady  preconditioning  for  the  blended  and 
AUSM+i/p'  schemes.  The  vorticity  profiles  and  inner  iteration  convergence  are  plotted  in  Figure  4.  At  this  much 
smaller  time  step,  the  steady  preconditioning  case  and  AUSM+i/p'  show  essentially  flat  convergence  while  the 
baseline  unsteady  preconditioning  and  blended  scheme  show  rapid  convergence.  From  Figure  4(a),  it  is  observed 
that  the  accuracy  for  the  baseline  unsteady  preconditioning  scheme  deteriorates  very  dramatically  while  the  blended 
scheme  continues  to  provide  an  accurate  vorticity  profile.  Despite  the  very  poor  convergence,  the  AUSM+wp* 
scheme  and  the  baseline  scheme  with  steady  preconditioning  give  as  accurate  a  profile  as  the  blended  flux 
formulation. 
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(b)  Inner  iteration  convergence 


Figure  4.  Vortex  solution  accuracy  and  convergence  for  a  time-step  of  CFLU  =  0.001. 


Thus  far  the  predicted  vorticity  field  was  considered  as  the  metric  for  evaluating  the  solution  accuracy. 
However,  while  the  pressure  variations  may  be  very  small,  its  gradient  should  correspond  to  the  swirl  velocity  in  the 
vortex.  Therefore,  it  is  instructive  to  evaluate  the  resulting  pressure  contours  for  the  difference  flux  procedures.  The 
vorticity  and  pressure  contours  for  the  various  schemes  are  plotted  in  Figure  5  and  Figure  6,  respectively.  The  flux 
difference  scheme  with  steady  preconditioning  gives  accurate  results  for  the  vorticity  field  as  seen  in  Figure  5(a), 
however,  the  pressure  field  in  Figure  6(a)  is  very  inaccurate.  This  is  not  unexpected  since  the  steady 
preconditioning  gives  the  correct  scale  for  the  velocity  equation  but  provides  an  incorrect  scale  for  the  acoustic 
eigenvalues.  The  flux  difference  scheme  with  unsteady  preconditioning  gives  poor  solutions  for  both  the  vorticity 
and  pressure  fields  in  Figure  5(b)  and  Figure  6(b)  due  to  large  dissipation  in  the  velocity  equation. 


(a)  Flux  Difference  with  (b)  Flux  Difference  with  (c)  Blended  Scheme  (d)  AUSM+i/p  ’  Scheme 

Steady  Preconditioning  Unsteady  Preconditioning 

Figure  5:  Vorticity  field  for  various  flux  schemes  for  a  time  step  of  CFLU=0.001. 


The  blended  flux  difference  scheme  gives  the  correct  solution  for  both  the  vorticity  and  pressure  field  as  shown 
in  Figure  5(c)  and  Figure  6(c),  since  they  use  different  preconditioning  scales  for  the  two  variables:  unsteady  for  the 
pressure  wave  and  steady  for  the  velocity  field.  However,  there  appears  to  be  one  notable  discrepancy  in  that  the 
pressure  field  is  not  circular  to  match  the  vortex  but  is  instead  distorted  to  be  a  rhombus.  The  cause  of  this  distortion 
is  not  fully  understood  but  may  be  due  to  the  difference  in  magnitude  for  the  cross-dissipation  terms  involving  the 
velocity  components  and  the  pressure  field  (e.g.,  AuAp  versus  AvA p). 

The  vorticity  and  pressure  fields  computed  using  the  modified  AUSM+w/U  scheme  can  be  seen  in  Figure  5(d) 
and  Figure  6(d).  Accurate  solutions  are  realized  for  both  the  velocity  and  pressure  fields.  In  particular,  the  pressure 
field  does  not  show  the  distortion  that  was  evident  from  the  blended  scheme  results.  Moreover,  the  circular  contours 
in  the  pressure  field  accurately  match  the  contours  of  the  velocity  field.  The  blended  flux-difference  and  the 
AUSM+i/p '  schemes  both  show  odd-even  oscillations  in  the  pressure  field.  The  source  of  this  instability  is  not  clear 
but  may  be  related  to  the  inconsistency  between  the  schemes  used  for  the  flux  for  the  implicit  and  explicit  sides. 
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(a)  Flux  Difference  with 
Steady  Preconditioning 


(b)  Flux  Difference  with 
Unsteady  Preconditioning 


(c)  Blended  Scheme 


(d)  AUSM+i/p '  Scheme 


Figure  6:  Gauge  pressure  field  for  various  flux  schemes  for  a  time  step  of  CFLU=.001. 


In  summary,  both  the  blended  flux  difference  and  the  AUSM+i/p  ’  schemes  are  able  to  provide  accurate  solutions 
of  the  vorticity  field  and  are  insensitive  to  the  time  step.  Both  schemes  are  a  substantial  improvement  over  previous 
baseline  flux  difference  procedures  that  showed  inconsistency  between  convergence  and  accuracy  with  the  solution 
accuracy  in  particular  deteriorating  rapidly  with  unsteady  preconditioning.  An  evaluation  of  the  pressure  contours 
revealed  that  both  the  blended  flux  difference  and  the  AUSM +up  ’  schemes  were  able  to  provide  the  correct  pressure 
depression  in  the  vortex.  However,  the  solutions  for  the  blended  flux  difference  scheme  show  a  distortion  of  the 
pressure  field  shape.  In  contrast,  the  pressure  field  computed  using  the  AUSM+w/f  scheme  maintained  the  circular 
pressure  field.  The  only  drawback  of  the  AUSM+i//?  ’  formulation  appears  to  be  less  than  ideal  convergence  at  small 
time  steps  and  this  is  probably  related  to  the  inconsistent  discretizations  on  the  LHS  and  RHS.  This  is  potentially  an 
area  requiring  improvement  in  terms  of  efficiency  of  the  procedure. 


C.  Simulations  of  Flow  in  a  Tube  with  Oscillating  Back  Pressure 

The  next  test  problem  investigated  was  inviscid  flow  in  a  1-D  tube  where  the  back  pressure  was  fluctuated  and 
the  inlet  pressure  was  held  constant  as  given  by 

Px=L=Pb(l+esmwt)  (31) 

P?=0  =  constant 


The  analytical  solution  for  this  scenario  is  given  below: 
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The  velocity  field  is  uniform  in  space  but  fluctuates  in  time,  while  the  pressure  field  has  a  linear  slope  spatially  from 
the  inlet  to  the  exit  while  also  fluctuating  in  time.  The  key  parameter  of  interest  here  is  the  Strouhal  number  Q;  as 
Q  increases  the  unsteady  time  scales  become  dominant  and  the  baseline  flux  difference  procedures  provide 
inaccurate  solutions  that  vary  with  the  preconditioning  parameter.  This  problem,  therefore,  would  be  a  good  test  for 
the  validity  of  both  the  blended  flux  difference  and  the  modified  AUSM+wp  ’  scheme.  The  objective  is  to  be  able  to 
get  an  accurate  solution  independent  of  the  preconditioning  parameter  and  frequency. 

In  the  present  calculations,  the  mean  Mach  number  in  the  tube  was  specified  to  be  0.001  and  a  5  percent 
fluctuation  amplitude  was  applied  at  the  exit  with  100  inner  iterations  for  the  dual  time  iteration.  The  frequency  and 
the  time  step  were  varied  over  a  broad  range  to  test  the  robustness  and  accuracy  of  the  flux  procedures  over  a  broad 
range  of  conditions.  The  steady  preconditioning  Mach  number  cut-off  was  specified  as  0.001  while  the  unsteady 
Mach  number  scale  was  set  to  1.0.  The  time  step  was  specified  as  2  n/(PPW Q)  where  PPW  denotes  the  points  per 
time  period. 

For  the  benign  scenario  where  the  fluctuation  frequency  is  low,  Q  =  1,  and  the  time-step  is  large  with  10  points 
per  time  period,  the  baseline  flux  difference  scheme  with  steady  and  unsteady  preconditioning,  the  blended  flux 
difference  scheme,  and  the  AUSM+up’  scheme  all  produced  pressure  and  velocity  profiles  that  matched  the  exact 
solution.  Essentially  in  this  case  the  flow  behaves  like  a  quasi- steady  problem  with  the  time  variation  being  a  series 
of  steady  variations.  The  convergence  history  shows  very  little  difference  between  steady  and  unsteady 
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preconditioning  for  the  various  schemes  and  is  consistent  with  this  quasi -steady  picture. 

The  results  for  the  more  difficult  case  in  which  the  fluctuation  frequency  is  high,  Cl  =  100,  and  the  time-step  is 
small  with  1000  points  per  time  period  are  presented  in  Figure  7.  For  the  flux  difference  scheme  with  steady 
preconditioning,  the  pressure  response  becomes  very  unstable  and  the  amplitude  overshoots  the  exact  solution  for 
velocity  and  pressure  as  seen  in  Figure  7(a)  and  Figure  7(b),  respectively.  Moreover,  an  excessive  phase  error  is 
also  evident.  The  pressure  residual  for  the  steady  preconditioning  scheme  is  shown  in  Figure  7(d)  and  shows  no 
convergence  during  the  inner  iterations.  This  is  consistent  with  the  incorrect  flow  results.  The  solutions  computed 
using  the  blended  flux  difference  scheme,  and  the  AUSM+up’  gave  essentially  the  same  results  that  matched  the 
exact  solution  very  well.  Therefore,  only  the  solution  profiles  predicted  by  the  modified  AUSM+up’  scheme  are 
shown  in  Figure  7.  The  velocity  profile  matches  the  exact  solution  perfectly.  However,  an  initial  fluctuation  around 
the  exact  solution  can  be  seen  for  the  pressure  transient  in  Figure  7(b).  This  fluctuation  decays  very  quickly  and  the 
source  of  the  error  is  not  clear  at  this  point.  The  convergence  history  plotted  in  Figure  7(c)  and  Figure  7(d)  for  the 
velocity  and  pressure  residuals  show  that  the  AUSM+up’  show  good  convergence  slopes  in  general  for  both  the 
velocity  and  pressure  residual.  In  contrast  the  steady  preconditioning  shows  very  poor  convergence  consistent  with 
the  fact  that  steady  preconditioning  is  inappropriate  for  high  frequency  pressure  fluctuations.  These  results  confirm 
that  for  low  Mach  number  acoustic  problems  that  the  dissipation  for  the  continuity  equation  (pressure  variable)  for 
both  the  AUSM+up’  and  the  blended  flux  difference  form  must  be  consistent  with  the  unsteady  preconditioning 
form. 


(c)  Velocity  residual  (d)  Pressure  residual 

Figure  7:  Velocity  and  pressure  history  and  residuals  for  a  high  oscillation  frequency  (Q  =  100)  computed 
with  a  small  time  step  (1000  PPW) 

D.  Simulations  of  Shock  Tube  with  Small  Pressure  Difference 

The  final  test  case  is  a  pure  acoustic  problem  with  low  velocities.  Here  we  model  a  shock  tube  with  a  very  small 
pressure  difference  across  the  diaphragm  that  generates  a  low  Mach  number  flow.  In  this  particular  case,  a  tube  of 
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length  1  m  was  chosen  and  the  diaphragm  was  placed  at  x  =  0.5.  The  pressure  on  the  left  side  of  the  diaphragm  is 
initialized  to  100028.04  Pa  and  the  pressure  on  the  right  side  is  initialized  to  100000  Pa;  therefore  the  pressure 
difference  is  28.04  Pa.  The  temperature  is  initialized  to  300  K  on  both  sides  and  the  initial  velocity  is  zero.  A  very 
weak  acoustic  waves  results  and  generates  a  low  Mach  number  flow  of  the  order  of  0.0001  in  the  contact  interface. 

In  the  calculations  presented  here,  the  time  step  was  kept  constant  at  1 0  ps  and  the  calculations  were  done  for  a 
total  time  of  750  ps.  One-hundred  inner  iterations  were  used  at  each  time-step  and  the  unsteady  Mach  number  for 
unsteady  preconditioning  was  set  to  1.0.  The  results  of  the  previous  problem  confirmed  that  dissipation 
corresponding  to  unsteady  preconditioning  must  be  employed  when  simulating  low  Mach  number  acoustic 
problems.  Therefore,  the  flux  difference  procedure  with  steady  preconditioning  will  not  be  considered  here. 

Figure  8  shows  the  solution  predicted  using  the  baseline  flux  difference  scheme  with  unsteady  preconditioning 
and  the  blended  flux  difference  procedure.  Recall  that  the  blended  scheme  uses  unsteady  preconditioning  for  the 
pressure  waves  and  steady  preconditioning  for  the  velocity  and  temperature  equation.  For  the  solution  profiles 
shown  in  Figure  8,  the  steady  preconditioning  parameter  is  set  to  Mi  =  0.01  and  the  unsteady  Mach  number  scale 
was  set  to  Mu  =  1.0.  It  can  be  seen  in  Figure  8(a)  and  Figure  8(b)  that  the  acoustic  wave  and  the  contact  interface 
have  propagated  to  the  correction  location  and  the  inner  iteration  convergence  history  shown  in  Figure  8(c)  appears 
to  be  good  for  both  flux  difference  procedures.  However,  closer  inspection  shows  that  the  blended  flux  difference 
procedure  produces  small  fluctuations  in  the  pressure  field  at  the  contact  interface  which  translates  to  substantial 
fluctuations  in  the  velocity.  The  level  of  the  fluctuations  is  sensitive  to  the  disparity  between  the  steady  and 
unsteady  preconditioning  in  the  blending  formula;  the  oscillation  amplitude  increases  as  the  steady  cut-off  Mach 
number  is  reduced  and  the  oscillation  amplitude  decreases  as  the  steady  cut-off  Mach  number  is  raised.  If  Mi  is 
increased  to  1.0  then  the  blended  flux  procedure  is  identical  to  the  flux  difference  scheme  with  unsteady 
preconditioning  which  does  not  produce  any  oscillations. 


(b)  Mach  number  profile  (c)  Inner  iteration  convergence 

Figure  8:  Shock  tube  solutions  for  the  blended  flux  difference  schemes  with  Mu=1.0  and  Ms=0.01 


The  profiles  predicted  using  the  AUSM+i/p  schemes  are  shown  in  Figure  9  for  two  different  formulations.  The 
profiles  given  by  the  blue-line  include  the  modification  to  the  mass  flux  dissipation  term  as  given  by  Eqn.  (27), 
however,  does  not  include  the  modification  to  the  pressure  dissipation  term  as  given  by  Eqn.  (29).  It  can  be 
observed  from  Figure  9  that  this  formulation  produces  fluctuations  near  the  contact  surface  similar  to  what  was  seen 
from  the  blended  flux  difference  procedure.  Additional  calculations  were  done  with  larger  pressure  ratios  that 
increase  the  flow  Mach  number  and  it  was  found  that  the  oscillations  dissipate  when  the  pressure  difference  is  large 
enough  to  generate  a  Mach  number  of  the  order  of  0.1  at  the  contact  interface.  This  provides  conclusive  evidence 
that  this  numerical  instability  is  related  to  acoustic  propagation  at  low  Mach  numbers.  In  particular,  these 
oscillations  appear  to  be  a  fundamental  manifestation  of  inadequate  dissipation  in  the  pressure  flux  that  arises  for 
low  Mach  number  flows  where  the  velocity  arises  from  the  propagating  pressure  pulse  and  appear  for  both  blended 
flux  difference  and  AUSM  family  of  schemes. 

The  derivation  of  the  additional  pressure  dissipation  in  Eqn.  (29)  was  in  fact  motivated  by  this  test  case  which 
pointed  out  the  inadequacy  of  the  pressure  dissipation  term  in  both  flux  differenced  and  AUSM  schemes.  The 
velocity  dissipation  term  in  the  pressure  flux  for  the  AUSM+i/p  formulation  goes  to  zero  as  the  Mach  number 
becomes  low.  However  this  dissipation  may  be  required  for  reducing  the  pressure  fluctuations  for  low  Mach 
number  acoustic  waves  and,  therefore,  the  unsteady  Mach  number  parameter  was  included  in  Eqn.  (29)  for  the 
definition  of  fa.  In  this  so-called  AUSM+w  ’p  ’  scheme,  the  addition  of  the  unsteady  Mach  number  scale  in  the 
velocity  dissipation  term  means  that  it  retains  a  substantial  value  for  unsteady  low  Mach  number  flows.  The  effect 
of  this  dissipation  is  illustrated  in  Figure  9  as  indicated  by  the  profiles  given  by  the  green  line.  The  addition  of  the 
velocity  dissipation  allows  for  solutions  that  are  smooth  in  both  the  pressure  and  Mach  number  profiles. 
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Furthermore,  the  convergence  is  slightly  improved  and,  therefore,  provides  a  scheme  that  works  in  an  accurate  and 
efficient  manner. 


Figure  9:  Shock  tube  solutions  for  the  AUSM+up  schemes  with  Mu=1.0 

In  summary,  simulations  for  low  Mach  number  shock  tube  configurations  which  represent  a  pure  acoustic  wave 
propagation  that  generates  a  very  low  Mach  number  flow  indicates  that  both  blended  flux  difference  and  the  original 
AUSM+wp/SLAU  schemes  have  some  fundamental  deficiencies.  While  the  location  of  the  acoustic  front  was 
accurately  captured,  the  pressure  and  velocity  in  the  contact  interface  showed  fluctuations  that  were  a  function  of  the 
local  Mach  number.  For  very  weak  acoustic  waves  where  the  flow  Mach  number  was  0.0001  substantial  oscillations 
were  observed  which  dissipated  as  the  Mach  number  rose  to  0.1  (for  a  larger  initial  pressure  difference).  The 
fundamental  deficiency  in  the  dissipation  was  rectified  in  the  AUSM+i/p'  scheme  by  adding  a  pressure  dissipation 
that  is  enforced  only  for  unsteady  low  Mach  number  flows  and  represents  a  fundamental  advancement  to  the  AUSM 
family  of  schemes.  Unfortunately  for  blended  flux  difference  schemes  there  is  no  clear  remedy  apart  from  going  to 
a  pure  unsteady  preconditioning  procedure  (i.e.  no  blending)  which  obviously  is  not  acceptable  for  the  broader  class 
of  unsteady  low  Mach  number  flows  that  were  tested. 

IV.  Concluding  Remarks 

A  generalized  preconditioning  framework  that  is  both  accurate  and  efficient  for  unsteady  low  Mach  number 
flows  is  presented.  It  rectifies  the  deficiencies  of  standard  steady/unsteady  preconditioned  flux  difference  schemes 
for  simulating  unsteady  low  Mach  number  flows  and  has  been  shown  to  apply  to  a  broad  class  of  problem 
encompassing  both  hydrodynamic  and  acoustic  driven  unsteady  flows.  Two  classes  of  schemes  were  considered: 
flux  difference  schemes  with  blended  steady  and  unsteady  preconditioning  and  AUSM  family  of  flux-splitting 
schemes.  For  flux  differenced  schemes,  generalized  “blending”  methodologies  were  developed  (extending  earlier 
work  by  Potsdam  et  al.1)  wherein  “unsteady”  preconditioning  is  used  for  the  pressure  wave  propagation,  while 
“steady”  preconditioning  is  used  for  the  convected  scalars  that  propagate  at  the  fluid  velocity.  Two  well-known 
AUSM  schemes  were  analyzed  initially;  the  AUSM+i/p  by  Liou5  and  the  SLAU  scheme  by  Shima  and  Kitamura7,8. 
The  SLAU  scheme  was  shown  to  be  equivalent  to  the  standard  AUSM+  formulation  (i.e.,  without  the  additional 
velocity  and  pressure  dissipation  proposed  for  low  Mach  number  flows).  A  more  generalized  formulation,  the 
AUSM+i/p'  scheme,  that  provides  a  unified  framework  for  unsteady  and  steady  flows  using  an  unsteady  Mach 
number  parameter  was  developed;  the  unsteady  Mach  number  parameter  which  can  be  specified  as  a  local  function 
of  the  flow  physics  is  shown  to  provide  the  right  selection  of  dissipation  in  the  mass  flux  term.  Furthermore  a  new 
velocity  dissipation  term  was  developed  in  the  AUSM+i/p'  formulation  for  acoustic  flows  to  suppress  spurious 
oscillations  at  very  low  Mach  numbers. 

Both  the  blended  flux  difference  schemes  and  the  modified  AUSM+wp'  schemes  were  tested  over  a  wide  range 
of  unsteady  test  cases  that  encompass  both  hydrodynamic  and  acoustic  unsteadiness.  Both  schemes  gave  superior 
results  with  both  the  pressure  wave  and  velocity  propagation  being  captured  accurately  independent  of  time  step  and 
Strouhal  number  while  providing  adequate  convergence  for  the  inner  iteration.  However  there  were  two  notable 
exceptions:  1)  for  multi-dimensional  flows  such  as  the  vortex  propagation  problem  the  blended  flux  difference 
schemes  distorted  the  pressure  field  while  the  AUSM+i/p'  scheme  did  not,  and  2)  for  small  pressure  jump  shock 
tube  problem,  oscillations  in  the  contact  surface  were  suppressed  in  the  AUSM+i/p'  scheme  by  the  addition  of  a 
velocity  dissipation  term  while  no  remedy  was  available  for  the  blended  flux  difference  procedure.  Therefore,  the 
modified  AUSM+wp  ’  and  AUSM +u  p  '  schemes  are  considered  the  superior  scheme  for  unsteady  low  Mach  number 
flows.  We  note  here  that  the  additional  velocity  dissipation  term  in  AUSM+w  p  ’  is  needed  only  for  suppressing  the 
oscillations  in  the  shock-tube  problem.  In  all  other  cases,  the  modified  AUSM+ipf  scheme  proved  to  be  the  best 
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choice.  Resolving  this  discrepancy  will  be  addressed  in  future  work.  Additionally,  improvement  of  the  inner 
iteration  convergence  for  the  AUSM+i/p'  scheme  also  deserves  further  scrutiny.  Here,  the  convergence  may 
currently  be  hindered  by  the  inconsistency  between  the  LHS  (implicit  matrix)  and  RHS.  Finally,  future  work  will 
also  focus  on  extending  the  methodology  to  frameworks  that  have  both  rotational  and  inertial  frames  as  well  as  on 
developing  generalizing  the  unsteady  preconditioning  parameter  to  consider  local  and  global  Strouhal  numbers. 
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